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Abstract. The thermally-driven formation and evolution of vertex domains is 
studied for square artificial spin ice. A self consistent mean field theory is used to show 
how domains of ground state ordering form spontaneously, and how these evolve in the 
presence of disorder. The role of fiuctuations is studied, using Monte Carlo simulations 
and analytical modelling. Domain wall dynamics are shown to be driven by a biasing 
of random fluctuations towards processes that shrink closed domains, and fluctuations 
within domains are shown to generate isolated small excitations, which may stabilise as 
the effective temperature is lowered. Domain dynamics and fluctuations are determined 
by interaction strengths, which are controlled by inter-element spacing. The role of 
interaction strength is studied via experiments and Monte Carlo simulations. Our 
mean field model is applicable to ferroelectric 'spin' ice, and we show that features 
similar to that of magnetic spin ice can be expected, but with different characteristic 
temperatures and rates. 
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1. Introduction 

Artificial spin ices [Ij are constructed as finite arrays of elongated magnetic dots 
whose magnetisations are assumed to be well approximated by Ising spins. The 
stray magnetostatic fields of each island mediate interactions, which are frustrated by 
geometry. These systems are called 'ices' because minimisation of the magnetostatic 
energy leads to behaviour resembling that governed by the ice rule for the ground state 
of solid water, i.e., two-in, two-out spin configurations [2j. 

The geometry of artificial spin ices can be controlled to a large extent [3H7], 
and resulting magnetic configurations can be imaged directly using techniques such 
as MFM [DIHHIlj, PEEM [SI [121 [13], and electron microscopy Most studies 

have been made at room temperature with relatively large, thermally stable magnetic 
elements. This thermal stability is enforced by choosing island volumes large enough 
that the barrier to magnetisation reversal, primarily associated with shape anisotropy, is 
much higher than room temperature, so that dynamics can only be induced by applying 
a magnetic field. Spatial studies are then made with the system in a steady state 
configuration. 

Because the focus to date has been on athermal systems, very little is known about 
how two dimensional spin ice arrays respond to thermal fluctuations. Recent reports 
provide evidence that there are magnetic features visible in as-grown arrays that form 
during early stages of growth [10], namely long range ground state ordering, with well 
defined domain walls and small clusters representing configurational excitations above 
the ground state. The origin of the ordering and excitations was argued to be thermal. In 
another recent work, thermal loss of macro-spin order was reported in a square artificial 
spin ice patterned (5-doped Pd(Fe) film [T5] . 

From a theoretical perspective, simulation studies [16] of an analogue of square ice 
constructed from charged colloids in double-well optical traps reveal a freezing transition 
as the temperature is reduced relative to the height of the barriers between pairs of 
wells. The same authors also find [17] that strong thermal noise is required to alter 
the hysteretic behaviour of the system. Recently, Monte Carlo simulations of magnetic 
square ice have been presented [T^, which indicate a sharp peak in the specific heat 
and a peak in the density of closed loops of spins flipped against the ground state at 
approximately the same temperature. 

In this paper we discuss configurational dynamics of square artificial spin ices in 
terms of macro-domain growth and boundary movements at finite temperatures. The 
picture we present is one in which domain boundaries can flow, and also form channels 
over which magnetically charged vertices move. We find that a domain of ground state 
configuration should, in the absence of disorder, spontaneously grow until it fills the 
array in a finite system. We show also that thermal fluctuations can accelerate domain 
boundary movement, and create clusters similar to those observed experimentally. 

We restrict our attention to square array geometries in this work. Spin 
configurations can be completely specified by vertex arrangements, and for the present 
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Figure 1. (a) Examples of each of the four vertex types. The symbols in the centres 
of vertices I-III are used later in the text to indicate the vertex type, (b) Geometry of 
the square artificial spin ice in this paper. Finite range interactions are assumed, and 
indicated on the diagram. 



work, we use spin or vertex descriptions interchangeably for configurational and energy 
states. The concept of vertices provides a useful nomenclature for discussing spin 
configurations [H [HI [161 [191 [20] , effective temperatures of field-driven demagnetisation 
[211122] and dynamics [m[23ti25] . We will use vertex descriptions almost exclusively 
in what follows. In the square lattice there are sixteen distinct vertex configurations, 
which can be classified into four groups. Type I — IV on the basis of magnetic charge 
and total moment. The groups are ordered on the basis of energy. Examples of each 
type are shown in figure [T](a). The ground state corresponds to a complete chessboard 
tiling of the two Type I vertices, and is two-fold degenerate. 

The organisation of the paper is as follows. In the next section we describe general 
features of domain dynamics in the presence of quenched disorder, using a mean field 
theory, in reference to conventional magnetic field and zero field cooling experiments. 
We then discuss effects of thermal fluctuations and strength of coupling using Monte 
Carlo simulations. These results are compared to experiments. We also discuss how 
thermal fluctuations can create small cluster excitations above the ground state, which 
can be understood using statistical arguments. Throughout this work, we see that Type 
I domain formation and growth are general phenomena, a result which is corroborated 
by simulation results for artificial 'spin' ice in ferroelectric media, which we discuss in 
an appendix. 

2. Magnetisation Processes in Mean Field Approximation 

Magnetic elements are approximated as block Ising spins, and arranged on a finite 
two dimensional square lattice, as shown in figure [H^b). The importance of long range 
dipolar interactions for correlations is a topic of investigation [T3l[T9l|26ll2Z], but our 
simulations suggest that shortcomings of the finite range interaction approximations 
are significant only in idealised perfect systems. In the models used throughout the 
main text of this paper, we consider only the three nearest neighbour interactions: Jn, 
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Jnn, and Jnnn in a poiiit dipole model. The square array geometry can be thought of as 
two sets of parallel lines of elements, with the sets aligned at 90° to each other. The Jnn 
accounts for interactions along each hne, and the Jnnn couples elements from adjacent 
parallel lines. The Jn couples elements from the two different sets. These couplings are 
indicated schematically by dashed lines in figure [^b). In the appendix, we give results 
for a ferroelectric ice where all point dipoles in the array are summed over, and find 
qualitatively similar results to those presented in this section. 

The energy of a magnetic element is determined by the local interactions and any 
applied fields, H. Denoting the magnetic moment at site i as m^, reversal of mj will 
occur when 



Here, /is is the Bohr magneton, (mj) is the mean field thermal average moment at site 
i, and ec is an energy barrier to reversal. We suppose that the energy barrier to reversal 
is proportional to the anisotropy responsible for the Ising like behaviour of the magnetic 
elements. In the mean field approximation this is 



where K represents an anisotropy barrier. Each magnetic element is assumed to behave 
as a single 'soft' macrospin, with a total moment that depends upon temperature. 
We approximate the temperature behaviour of an element's moment by the Langevin 
function L{x): 



f3 is the inverse temperature 1/(A;bT) where ks is the Boltzmann constant. 

Previous studies have shown that, in field-driven dynamics, a perfect system allows 
nucleation of Type I domains on a Type II background through Type III vertex dynamics 
that start at array edges where moments have fewer neighbours and therefore lower 
reversal thresholds [21]. Local variations in the coupling or reversal barrier parameters 
can facilitate nucleation processes to start inside the array, leading to less sensitivity to 
the array boundaries [20l|28l|29]. In a simple model for switching barriers, the barrier 
height is proportional to the element volume. One would expect that during early 
stages of deposition of material, the element thicknesses are distributed over a range, 
leading to a spread in switching field values. Furthermore, significant distributions in 
switching fields have been observed experimentally in fully-grown, athermal artificial 
spin ice p^[29ti32] and switching field disorder has been shown to give similar outcomes 
in numerical simulations as other disorder types [33j. We use this as motivation for 
describing disorder in terms of the reversal barriers. We assume the barriers vary 
randomly in an inverval A centred about K, that is, [K — A/2, K + A/2]. 

The mean field model allows one to examine the stability of the ground state relative 
to temperature but cannot capture correlations such as those involved in avalanche 
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processes, and such processes will be discussed later. Initially, all moments are assumed 
to be of unit magnitude in some specified configuration. In the present model, an 
element i is picked at random, and the time averaged reduction of its magnetisation is 
calculated according to ([3]). If the magnetisation state becomes unstable, that is, if ([1]) 
is satisfied, then rrii is set to — mj. The process is repeated by choosing a new moment 
at random, and the iterations continue until a steady state is found such that no rrii 
value changes by more than 10~^. 

Reduced units used throughout the paper are defined so that |mj| < 1, Jn = 1-5 Jo, 
Jnn = 0.7Jo, and J„„„ = O.SJq. The reversal barrier has been chosen to be, in reduced 
units, K = lOJo- Arrays consist of 640 vertices, and array boundaries are taken to be 
'open' boundaries for which the edge elements have only three nearest neighbors, and 
corner elements have two nearest neighbours, as shown in figure [TJ^b). 

The stability of the Type I ground state can be studied by starting the system at 
low temperature with a complete Type I tiling, and then increasing the temperature. 
The critical temperature can be calculated within mean field theory from the average 
magnetisation of an element in a Type I tiling, which obeys 



where /iioc = 2(2Jn — Jnn + Jnnn) + . The critical temperature Tc, is given by 
= /iioc/3, to first order approximation. Using the parameters listed above, 
Tc = 3.4To, where To = Jo/ks- 

Because changes in magnetic configurations are driven by instabilities only, in the 
mean field model, the transition is very sharp with a well defined critical temperature in a 
perfect system without disorder. Disorder broadens the transition region in temperature, 
and leads to formation of domains of the two different possible Type I tilings separated 
by Type II and III chains. Examples are shown in figure [2] for Type I populations ni as 
functions of temperature for different disorders A. The population shown is the sum of 
both Type I tiling possibilities. Other authors also report a broad melting transition in 
a thermal ice prepared in a magnetised (Type II vertex tiling) state [15] . 

With disorder, a sharp initial reduction occurs in ni near the critical temperature, 
corresponding to the nucleation and growth of domains. This is followed by a long 
tail in which the total population of Type I vertices approaches the value ni = 1/8, 
corresponding to the high temperature limit in which each of the 16 possible vertices 
appear with equal probability. Note that the transition temperature goes to the expected 
value of 3.4To as A — > 0. 

Also shown in figure [2] is the Type I population that evolves from a random 
vertex configuration starting at low temperature. This is analogous to the 'zero field 
cooled' state of a ferromagnet that, initially at some temperature above the ordering 
temperature, has been quenched to low temperature. In the mean field model, a state 
with large Type I domains immediately evolves at low temperature if K is sufficiently 
small relative to the J's. The rii populations grow with increasing temperature as other 
vertices become unstable. 
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Figure 2. Thermal reduction of Type I population starting from a random 
configuration (solid lines) and starting from a perfect ground state (dashed lines). 
Disorder is introduced through random variations in individual barriers to reversal. 
Thermal effects are calculated in a self consistent mean field approximation. The 
total population of all Type I vertices is shown as a function of reduced temperature 
for different amounts of disorder. For each disorder value, the evolution is traced 
starting from low temperature for a uniform single domain of Type I, and a random 
configuration. The transition temperature approaches the expected value S.GTq in the 
limit of no disorder. Disorder broadens the transition considerably. 

This growth of domains occurs through the ehmination of small domains and 
clusters. An example is shown in figure El where vertex configurations at two steps 
during the numerical iteration are shown. The temperature is To, and disorder A = 0.1. 
Green squares represent Type I vertices. The blue arrows are Type II vertices, and their 
direction indicate the orientation of the local Type II moment. The light and dark red 
crosses represent the two charge flavors of the Type III vertices. 

The first image (left) is taken after 100 iterations after an initially random 
configuration of moments. Already a well defined Type I domain structure has formed. 
The domains are separated by chains composed of Type II and III vertices with net 
magnetic moment and charge, and the energy of their formation acts as a surface tension 
on the Type I domains. For this reason, the chains tend towards straight in the limit 
of no disorder, since corners lengthen the chain and can involve Type Ill's. The second 
image (right) shows the evolution after 100 additional iterations. The longest chains have 
shifted slightly, and positions of the Type III vertices within the chains have changed. 
Most notably, the smallest domains have vanished, and other small domains are reduced 
in size. Note that Type IV vertices are highly unstable, and do not survive past the 
first few iterations. 

The reduction in size and eventual annihilation of small domains is the result 
of an energetic biasing of random domain wall motion towards motion that makes 
walls shorter by moving closed walls inwards. Dynamics on a Type I/Type II vertex 
background - such as domain wall motion - can be described in terms of Type III vertex 
motion pT | [23 | l2l]. There are therefore two energy costs associated with dynamics: the 
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Figure 3. Example configurations during numerical iteration of the self consistent 
mean field algorithm. Green regions are Type I domains. Blue arrows indicate Type 

II vertices, and point in the direction of the vertex net moment. Red crosses are Type 

III vertices. The left panel shows a configuration found after 100 iteration loops begun 
from a random configuration at temperature To and no disorder. The right panel is 
generated 100 iteration loops later. The box in each panel indicates a closed domain 
which has shrunk over the course of the evolution. 



cost of creating Type III vertices, and the cost of their propagation. This latter cost 
is approximately zero in processes where the vertex populations are conserved. For 
example, in figure Ht^a), Type III vertices on the domain wall can move at approximately 
zero cost by flipping the circled spins, to yield the configuration shown in figure |l][b) . 
Accordingly, the lowest-energy mode of domain wall motion is the propagation of Type 
III vertices that have formed at random positions during wall creation. These random 
fluctuations do not favour either growth or shrinking of domains on average. The 
existence of random fluctuations lead us to speculate that an analogy to creep motion [31] 
may be possible. 

The cost of Type III vertex creation depends on the local spin configuration, and 
is lowest at corners in the domain walls. For example, in figure ID^c), the easiest spin 
to flip is the circled spin at the domain wall corner, which can flip to create a Type III 
vertex pair. The resulting Type III vertices can propagate at approximately zero cost 
by moving along the domain wall, flipping a diagonal chain of spins on the inside of the 
wall. This process moves the domain wall inwards, as illustrated in figure ID^d). Such 
sequences of Type III creation and propagation are the domain wall motion mechanism 
with second-lowest energy cost. Thus, over time, the wall motion is biased so that 
despite random fluctuations, all closed domains will disappear. 

As a final comment on the mean field results, the local fields along the element 
magnetisations are weakened within the chains, compared to the fields in regions of 
lower-energy Type I vertices. In the mean field model this leads to a suppression of 
the magnetisation in the elements participating in the chains. The magnitude of the 
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Figure 4. (a, b) Random domain wall fluctuations, driven by Type III vertex motion. 
The circled spins in (a) can flip at approximately zero cost, to yield the configuration 
shown in (b). (c, d) Domain shrinking by Type III vertex creation and propagation. 
The corner of the domain wall provides a Type III nucleation site, shown as the circled 
spin in (c). Once nucleated, the Type III vertex can propagate at approximately zero 
energy cost, shrinking the darker blue domain, as seen in (d). The arrows representing 
spins are colour-coded light green and dark blue according to which Type I domain 
they belong to. Vertex types are indicated by symbols: Type I, II and III vertices 
represented as circles, arrows and crosses, respectively. 



suppression depends upon the temperature, and can be very large as one nears Tc. This 
effect is shown by the vector magnitude map in figure [5l The magnitude of element 
magnetisation is shown in grayscale, with black being unity and white being small. The 
domain pattern corresponds to the configuration shown in the second panel of figure [3l 
One sees that the moment is considerably reduced within the chains, and especially for 
small clusters. We note that instabilities occur when the moment of individual elements 
vanish, and these are most likely to occur in the Type II walls. 



3. Fluctuations and Disorder 



The energy differences between configurations are often small, even for configurations 
that are very different. For this reason, single spin flips can drive avalanches that lead to 
large configurational changes and strongly modify vertex populations as domains evolve. 
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Figure 5. The reduction of magnetic moment at finite temperatures is largest within 
Type II and III chains bounding Type I domains. Here a greyscale plot of the moment 
is shown for the second example configuration of figure |3l Dark signifies large moment, 
and white signifies small moment. The dark regions correspond to Type I domains, 
whereas the moment is dramatically reduced in the domain walls. 

To capture these dynamics, we turn to Monte Carlo simulations. 

The model is based as before on the point dipole approximation assumed for 
the mean field model. The form of the energy used is the same, but evolution is 
described using a heat bath algorithm, rather than a simple critical field defined by 
([3]). Additionally, we consider the case of a large Curie temperature so that the effective 
block moments are insensitive to temperature. In this limit the energy used for the heat 
bath algorithm is 

e = mj{^BH +^Kmj + ^Jj^kmk), (5) 

where now all magnetisations have fixed magnitude as in an Ising model, with \mi\ = 1. 
The barrier to reversal is determined by the anisotropy constant K. The parameters 
used are the same as those for the mean field model. The results shown below were made 
with 10, 000 Monte Carlo steps at each temperature and averaged over 20 realisations 
of disorder at each disorder strength. 

As in the mean field calculation, the role of disorder is largely to broaden the 
transition. Example heating curves are shown in figure [6], where thermal evolution of 
ground state and random tilings are shown, for different disorder strengths A. The most 
striking feature, in comparison to the analogous mean field results shown in figure [21 
is the smooth behaviour of the population at the transition, and the rounding in the 
transition region. 

The dependence on disorder is also less dramatic than in the mean field case. In both 
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Figure 6. Thermal fluctuation effects in an artificial square ice calculated using Monte 
Carlo, and shown as a function of temperature. The total population of all Type I 
vertices is shown as a function of reduced temperature for different amounts of disorder. 
For each disorder value, the evolution is traced starting from low temperature for a 
uniform single domain of Type I (dashed lines), and a random configuration (solid 
lines). At high temperatures the population tends towards 1/8, the value expected for 
a random sampling of the 16 possible vertex types. 




mean field and Monte Carlo simulations, the spatial distribution of switching barriers 
is weakly connected to the final distribution of domains and chains. An analysis of the 
distribution of K values shows that there is a propensity for chains to locate on strong 
pinning sites, where K is large, and for weak pinning sites to lie more frequently inside 
domains. This result can be expected because domains grow via thermally activated 
spin flips. Growth is accomplished through distortions of chains, and chains will pin at 
sites where K values are too large to allow thermally driven element magnetisation flips. 
Also, as a consequence, the final resulting domain configuration at a given temperature 
displays rough chains for large A values, and flat segmented chains for low A values. 

Differences between Monte Carlo and mean field results appear because of how 
thermal fluctuations drive domain growth. Fluctuations allow the system to explore 
dynamical pathways different to those dictated by a stability analysis, and are therefore 
quite significant in determining the final domain configuration. At temperatures below 
the transition, the dominant paths tend towards the lower energy, single domain state. 

An example evolution is shown in figure [TJ using the colour scheme defined in 
figure 131 Four configurations are shown at constant temperature T = To, starting with 
a random initial configuration. The first three configurations are taken after 1000, 2000, 
and 3000 steps. The fourth configuration (in the bottom right hand corner) is taken 
after 10, 000 Monte Carlo steps. The amount of disorder in this example was small 
relative to both K = IOJq and = l.bJo, with A = 0.1. One sees the slow growth of 
one Type I domain at the expense of another, through the motion and deformation of 
Type II and III chains. The dynamics is equivalent to that seen with the mean field 
model. Note that here also Type IV vertices do not survive the first few Monte Carlo 
steps unless interactions are weak (as with large spacings). 
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Figure 7. Different stages in evolution of a vertex configuration as calculated using 
Monte Carlo simulations. Green regions are Type I domains. Blue arrows indicate 
Type II vertices, and red crosses are Type III vertices. The first three panels show 
configurations after 1, 000, 2, 000, and 3, 000 steps for a temperature of Tq and disorder 
A = 0.1. The initial state was random. The final panel shows the configuration after 
10, 000 steps. Note that small clusters and domains disappear as large domains grow 
at the expense of smaller domains. 

In addition to our mean-field and Monte Carlo simulations, we have imaged vertex 
configurations in as-grown samples, and found domain structures similar to those in 
simulations. An example MFM image is shown in figure El where domain walls, which 
carry net moment are clearly visible on the approximately neutral background of Type 
I domains. The sample imaged here was previously studied in terms of ground state 
ordering and thermal excitations [10], and consists of an array of nominally 280 x 85 x 26 
nm^ Permalloy islands, with a 3 nm Ti buffer and a 2.5 nm Al cap. Full details of 
sample fabrication are given elsewhere [10] , but the key point is that islands are formed 
by deposition of Permalloy through gaps in a nanopatterned resist mask, so that island 
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Figure 8. MFM image of an as-grown square artificial spin ice, with lattice constant 
400 nm, relative to island dimensions of 280 nm x 85 nm. Red and blue colouring 
indicate magnetic poles, and domain walls and small excitations are clearly visible on 
a ground state background. Two examples of small cluster excitations are indicated 
by boxes. The image has been false-coloured using the software package WSxM [35] . 
Note that islands in this image are at 45° to the islands in other figures in this paper. 



thicknesses - and hence barriers to magnetisation reversal - grow over time. 

For the simulation, the values shown in the figure were obtained after a finite 
number of Monte Carlo steps, and the effect of temperature, relative to the local coercive 
field, is to control the rate at which ni changes. In the experiment, the deposition rate 
relative to temperature probably controls the rate at which ni changes. We suggest 
that domain configurations found in the simulation and the experiment can both bo 
interpreted as snapshots in time of a thermally driven evolution. We note that a 
difference between the simulations shown in this example and the experiment does exist, 
in that the as-grown samples support small clusters of flipped spins. We discuss this 
further in section [3l2l 

3.1. Spacing and interaction strength 

The preponderance of Type I vertices in domain structures is a consequence of the higher 
energies needed to form other vertex types, and the corresponding lower probability for 
their creation. The probability for reversing a moment depends upon the size of local 
flelds relative to thermal energy. This ratio also determines the rate at which domains 
can grow. 

An example from Monte Carlo simulations is shown in flgure[9](a). Here the Type 
I population is shown as function of interaction strength. The interaction strength 
is scaled by a factor s such that Jq, — )■ sJa/Jo for each of the three interactions a. 
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Figure 9. (a) The average number of Type I vertices as a function of interaction 
strength scaling s at a temperature of T/Tq = 1 for different strengths of disorder 
as simulated using Monte Carlo. Increasing disorder leads to greater stability for 
other vertex types relative to the ground state. The ni population tends toward 
unity as the simulation time is increased in all cases, but very slowly when disorder 
is strong, (b) Experimental results for Type I populations as determined from as- 
grown samples for arrays grown with different lattice spacings. The populations are 
shown as a function of where r is the lattice spacing, which is proportional to 

the interaction strength. The populations are the averaged over a number of images 
(typically 5 or 6) taken across each continuous pattern of 0.5 mm by 0.5 mm total 
area, and the errors are the standard errors over these. 



Temperature is fixed at To in all cases. Several different disorder strengths were studied 
and the simulations were run for 100, 000 steps at each spacing. The rii population tends 
to unity as the interaction strength increases, although the rate at which this occurs 
depends on disorder strength. For strong disorder, the rate is slowest. We note also 
that there were no significant Type IV populations observed, although the population 
did increase with increased disorder. 

We have also measured statistics of the vertex type populations for as-grown 
samples with different lattice constants. In these studies, the dimensions of individual 
elements were held constant at 270 x 115 x 25 nm^, but the lattice constant was varied 
from 400 nm to 600 nm. In this way the interaction strength should decrease, roughly 
like 1/r^ where r is the lattice constant. Results are shown in figure |9]|^b), where the ni 
population is given as a function of 1/r^. All samples in the series consist of 0.5 x 0.5 mm^ 
continuous arrays of Permalloy islands grown on a 2 nm Ti buffer, with no capping 
layer, and were fabricated in a single batch, to ensure that fabrication parameters were 
consistent across the series and so avoid problems with quantitative reproducibility (such 
as those mentioned in [TO]). 

Both figures |9]^a) and (b) show that the Type I populations increase with increasing 
interactions. In principle one should be able to extract estimates of interaction energies 
as a function of spacing from the experimentally measured ni. However, inspection of 
figure [U^a) shows that one also needs to know the amount of disorder before a prediction 
for ui can be made. The disorder can be expected to be strongly dependent on details 
of the sample growth and design [32] . 
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3.2. Fluctuations and clusters 

The role of disorder is subtle. On the one hand, disorder in K leads to randomness in 
the chains. On the other, disorder also facilitates fluctuations and leads to appearance 
of small clusters of Type II and III excitations, examples of which are illustrated in 
figure [TOl In Monte Carlo simulations, these fluctuations usually appear and disappear 
quite rapidly (in Monte Carlo step time), but can also lead to structures that may persist 
through several Monte Carlo steps. Experimentally, as-grown samples are seen to exhibit 
frozen-in fluctuations, as seen in figure [HI In previous studies, it was shown that the 
population of clusters decays approximately exponentially with their energy, which was 
given as evidence of thermally-driven processes occurring during sample growth [TO] . 
In simulations of field-annealed square ices realised in nanostructured superconductors, 
small structures can only exist independently of domain walls when sufficient disorder 
is present [20] . 

An analytical description of the nucleation, growth and decay of an isolated cluster 
gives insight into the distribution of experimentally observed clusters. As seen in 
figure Uni clusters can be thought of as connected lines of spin flips against the ground 
state. The lines may branch. In our analytical model, the growth and decay of clusters 
proceeds by single spin flips that extend or shorten these lines. For example, the left 
hand cluster shown in figure [TO] can evolve in to the right hand cluster, and vice- versa. 
For simplicity, we neglect disjoint clusters: large clusters cannot evolve into two smaller 
clusters. 

Under the assumption that clusters are excitations above the ground state that are 
nucleated by random thermal spin flips, we start with an initial perfect ground state 
and study the temperature-dependent evolution of small clusters of up to four spin flips. 
The clusters may be grouped by their topology (and hence energy). There are twenty 
five groups of such clusters: apart from the ground state, there is one distinct cluster 
of one flip, two clusters of two flips, five of three flips and sixteen of four flips. These 
groups include clusters that contain Type IV vertices, but these high-energy clusters are 
suppressed relative to those containing only Types II and III vertices. 

The processes of cluster growth and decay are transitions between the groups of 
clusters, and can be described by a master equation for the probability, P{A), of a 




Figure 10. Example 2 and 3 spin-flip cluster excitations on a Type I background. 
The arrows represent island moments, with the flipped spins shown as bold red arrows. 
The vertex configurations are also shown with Type I, II and III vertices represented 
as circles, arrows and crosses, respectively. 
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Figure 11. Population vs energy above the ground state (in units of u = /xqto^/ (Airr^), 
where m is the island moment and r is the lattice spacing) for cluster excitations of up 
to 4 spin flips. As indicated in the legend, data are shown from experiment (these were 
previously published in [10.), the master equations ^ and Boltzmann factors. The 
labels indicate the mnemonics used in [TU]. The labels 4+ and 4t are in parentheses 
to indicate that those clusters were not observed experimentally in |TD] . Error bars on 
the experimental data correspond to the square root of the population of each cluster 
type. 



cluster having topology A: 

dP{A,t)/dt = Xl(<^(^ ^ A)P{B,t) 

B 

- G{A B)u{A B)P{A, t)) , (6) 

where the sum runs over all cluster topologies that a cluster of type A can evolve to 
or from by a single spin flip. G{A — t- 5) is a multiplicity that takes into account the 
number of 'pathways' from A to B, and depends on the numbers of orientationally- 
distinct clusters with a given topology and the number of clusters of topology A that 
can evolve into a particular cluster of topology B. The rate z/(yl — ?■ 5) is given by the 
Arrhenius-law factor 

u{A ^B) = /exp(-(ee - edip)/A;BT). (7) 

Cc — Cdip gives the energy cost of a spin flip to grow or shrink a cluster on a Type I 
background, ec is an arbitrary anisotropy barrier whose exact value does not affect 
the steady state cluster probabilities. In this model, disorder in anisotropy barriers is 
neglected, and ec is constant always, edip is the energy of the spin prior to flipping, due 
to point-dipole interactions among near neighbours (Jn, Jnn and Jnnn)- The attempt 
frequency / does not affect the steady state probabilities. 

The twenty five coupled equations of the form ([6]) can be solved numerically, 
with an initial condition that the system is in the ground state: P{0,t = 0) = 1, 
P{i,t = 0) = 0,Vz 7^ 0, where represents the ground state. The steady state 
probabilities of all clusters that do not contain Type IV vertices, normalised to the 
experimentally observed population of the 1 excitation (that is, a single spin flip), are 
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plotted in figure [TT], as a function of their energy difference witli tlie ground state (as 
given in flO\). Tfie agreement between tlieory and experiment is remarkable. 

The figure shows the solution at fc^T = 2.96m, the temperature at which the 
exponential decay that best fits the theoretical results matches the exponential decay 
of experimental data. The energy scale u = iiom?' / (47rr^), where r is the lattice spacing 
and the moment m = MV, where M = 860 x 10^ Am^^ is the magnetisation and V is 
the island volume. The clusters that contain Type IV vertices are not shown, but apart 
from two exceptions they all have populations less than half of the lowest population 
shown. This is consistent with experiments [10], where such clusters are not seen and it 
was concluded such clusters were highly improbable. The exceptions are a Type III - 
Type IV - Type III line and a Type III - Type IV - Type II - Type III cluster, which 
both have predicted populations of ~ 2. 

For comparison, the figure also shows populations predicted from Boltzmann factors 
of the form g exp{—dE / ksT) . The degeneracy g of a cluster topology is determined by 
the number of distinguishable ways a particular shape can be rotated and reflected, as 
well as a factor of 2 (equal for all clusters) to account for the possibility of a global spin 
flip. dE is the cluster's energy above the ground state, taken from The ratio of 
interaction strength to temperature has been tuned to match the exponential decay to 
that of the experimental data, giving a ratio of ksT = 8.18m. 

In experiments, configurations are frozen in as the island volume becomes large 
enough that shape anisotropy barriers suppress island magnetisation reversal. The ratio 
of temperature to nearest-neighbour coupling can be compared for the master equations 
and the Boltzmann distribution by estimating the thickness at which freezing-in occurred 
from the two models. If we estimate the temperature during growth to be 350 K, then 
the nearest-neighbour interaction is 2.4 x 10^^^ J for the master equation and 5.9 x 10^^^ 
J for the Boltzmann factors. For Morgan et a/.'s 280 x 85 nm^ islands, the estimated 
interaction strengths correspond to thicknesses of 3.5 nm (master equation) and 1.7 nm 
(Boltzmann factors). Both these estimates are of the same order of magnitude as the 
estimate of ~ 1 nm given in [10]. 

Interestingly, while the Boltzmann factors and the solutions to ([6]) both agree well 
with the experimental data, the Boltzmann factors agree best for the excitations 1, 2L, 
3Z, 4Z while the solutions to (I6l) fit better for the other clusters. For example, the 
'pathway' 1 — )■ 2L — )■ 3U — )■ 40 leads to a stable configuration, 40. In ([6]), we have 
truncated the maximum cluster size to 4 flips, making 40 very stable. However, even 
without this approximation, the energy cost of adding an extra spin flip to 40 is high, 
and the cluster is stable. Thus, the 40 cluster is in some sense a 'sink' of probability, 
and has higher probability than would be expected from Boltzmann factors, which is 
captured by the master equations. 

Finally, we note that the master equations ([6]) explicitly assume an initial ground 
state configuration, and treats small excitations as having nucleated and grown on this 
background before being frozen by increasing island switching barriers. The clusters in 
the as-grown samples of Morgan et al. might also be interpreted as regions which have 
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fallen out of equilibrium during the effective cooling of an initially random configuration. 
While we have not explicitly tested such a possibility, the agreement of the master 
equations with experiment gives further evidence for the picture of nucleating and 
growing clusters. 

4. Conclusion 

Vertices in artificial spin ice are local configurations of magnetisations that can be 
thought of to some extent as classically emergent objects analogous to quasiparticles, 
able to move and interact according to well defined rules dictated by the lattice geometry. 
We have shown that domains of ground state vertices form spontaneously in the square 
lattice. Domain boundaries are defined by chains of Type II vertices, along which Type 
III vertices can appear, move and drive growth of one domain at the expense of another. 

We have shown that these domain dynamics occur with rates that depend on 
temperature and the strengths of interactions between elements. Effects of disorder have 
been studied, and shown to affect mainly the average size and growth rate of domains, 
while not modifying significantly the fundamental processes involved. A comparison 
to experimental results was obtained from as-grown samples for which interelement 
spacing was varied. The experiment and simulations are essentially different in that 
element thicknesses are changing with time in the experiment whereas temperature is 
changed in the simulation. Nevertheless, similarities between populations for the Type I 
ground states were observed between experiment and simulation, suggesting that the as- 
grown samples are behaving somewhat as frozen snapshots of magnetic ordering within 
the arrays occurring as described by the simulations presented here. 

In addition to domain wall motion, there is also the possibility that an element 
will flip within a domain. Reversal of a single element in a Type I configuration creates 
a pair of oppositely charged Type III vertices. Additional reversals lead to clusters of 
Type II and Type III vertices, that may persist for some time. We have calculated the 
probabilities of occurrence of small clusters as a function of their energy, and shown 
that one obtains a distribution in quantitative agreement with that reported recently 
from experiments [TO] . 

Lastly, we have shown that all features of domain formation and growth can be 
obtained with a generic mean field model. This leads us to suggest the possibility of 
creating artificial spin ice using ferroelectric media, and we have provided an example 
in the appendix using parameters appropriate to bulk PbTiOa. In this example we also 
showed that models using full dipole sums over a lattice of point dipoles produces the 
same qualitative results as a model using a severely truncated sum. 
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Appendix: Ferroelectric Ice 

Our discussion so far lias been within the context of magnetic spin ice, although our 
results have general applicability to a number of different systems. We illustrate this 
with a return to the mean field model, and discuss now domain growth within the 
context of Landau- Ginzburg theory for ferroelectrics. 

A second order Landau-Ginzburg model is used to calculate the electric dipole 
moment of an island i, given the effective electric field E at that point. Each dipole 
is assumed to lie along only one axis, parallel to the long axis of that island. The 
ferroelectric free energy density F at island position ri is given by: 

F{i) = '^P\^) + ^-P'{^)-^{i)■m, (8) 

where P is the electric polarisation, a = A(T — T^) and /5 are the Landau parameters for 
the material. P{i) is assumed to lie parallel with the element long axis. A point dipole 
approximation is again used with the electric dipole associated with island i denoted as 
Pj = VP{i), where V is the volume of an element and P(z) is assumed to be uniform 
across the element. 

Unlike the previous mean field theory, in this model we calculate the electric field at 
position i by superposing contributions from all other elements (approximated as point 
dipoles) in the square lattice. The number of elements used is sufficiently small that 
one can use the simple summation 



E(i') = V ^ f ^ {rj, - r,)[{ri, - Tj) ■ Pi] _ \ 



(9) 

Minimisation of ([H]) will, for a range of E values, give two solutions for pi with 
opposite sign. One can show simply from ([8]) that the energy barrier separating the two 
solutions disappears when 

> . 10 

Pi V27/31/2 ^ ^ 

This is the condition for stability of an element's polarisation orientation. When this 
condition is fulfilled, the polarisation is reversed. As in the previous mean field model, 
this procedure is iterated through a system of elements at a fixed temperature until a 
steady state is reached. 

Parameters used for the simulations are: Tc = 1100 K, A = 7.5 x 10^ C^^ m^ N K~^, 
and f3 = 2.4 x 10^ m^ N. These parameters give a polarisation in zero field and at 
room temperature P = \J\a/l3) = 0.5 C/m^, which is typical for PbTiOs in bulk [36] 
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Figure 12. The average size of a connected region (domain) of Type I vertices as a 
function of temperature as calculated using a self consistent Landau-Ginzburg model 
with parameters appropriate for PbTiOs. The error bars show one standard deviation. 



and is quite accurate for elements that are over 100 nm long. It should be noted 
that the Landau expansion (IHl) should strictly speaking be to sixth-order to most 
accurately model the ferroelectric phase transition [36] but we neglect the terms 
here for simplicity. Results are given for a 20 micron square sample with = 50 islands 
along one edge (4900 islands and 2401 vertices). A spacing of 400 nm between island 
centres is assumed. The starting configuration for an iteration is taken to be islands 
having polarisation with constant magnitude Pi = {^Y^'^, with random alignments along 
element axes. 

The energy barrier with these parameters for flipping an element polarisation is 
much larger than the thermal energy, for most temperatures below the Curie point, so 
reversal via thermal fluctuations occurs only in the vicinity of Tc. 

The general features observed in the magnetic system are found also for 
the ferroelectric system. Starting from a macroscopically upolarised state at low 
temperature, the effect of increasing temperature is to increase the number of Type 
I vertices. As before, domains of Type I's form, separated by chains of Type II and III 
vertices. Effects of disorder are also completely analogous. 

Finally, we note that domains grow at temperatures approaching Tc, as also found 
for the magnetic system. An example of average Type I domain size is shown in figure 
[T2] for three temperatures. The error bars correspond to the spread in sizes calculated 
over ten different initial configurations. 

It is useful to note that it is possible to define a relaxational dynamics based on 
the effective field — The iteration process can thus be related directly to a real 
time dynamics. In this sense, the numerical iteration method is a time evolution, and 
these simulations give a feeling for the effect that different cooling rates may have on a 
system. If the system is cooled faster than the islands' polarisation can respond, then 
a higher number of unfavourable Type II and Type III vertices will be frozen into the 
artificial ferroelectric ice. This is interesting because these domain walls form without 
the presence of disorder and have a density depending on the history of the system. 



